/*
 * Copyright 2010-2012 Susanta Tewari. <freecode4susant@users.sourceforge.net>
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

package genomemap.stat.tests;

import JSci.maths.statistics.ChiSqrDistribution;

import genomemap.data.Clone;
import genomemap.data.PhysicalMapService;
import genomemap.tagging.LinkingCloneStrategy;

import javautil.collections.ArrayUtil;
import javautil.io.LogUtil;
import javautil.lang.MathUtil;
import javautil.stat.desc.FrequencyUtil;
import javautil.stat.desc.SampleTotal;

import java.text.DecimalFormat;

/**
 * A strong inference framework was used to examine the distribution
 *  of repeats along chromosomes.
 *  clones are collected from the minimal tiling map i.e., probes in the order
 *  with their connecting clones.
 *
 *  H0:  The S,R,O clones are randomly distributed along a chromosome
 *  with constant probabilities of occurrence along the chromosome.
 *  H1:  The S,R,O clones constitute a stationary Markov Chain.
 *  A Chi-squared Test of H0 vs. H1 was performed (Basawa, I.V. and
 *  B.L.S. P. Rao (1980), Statistical Inference for Stochastic Processes,
 *  Academic Press, p. 62).
 *
 * @author stewari1
 */
public class RepDNATests {

    /** Field description */
    private final int m = 3;

    /** Field description */
    private final boolean debug = true;

    /** Field description */
    private final int[] SSCount;

    /** Field description */
    private final int[] SRCount;

    /** Field description */
    private final int[] SACount;

    /** Field description */
    private final int[] RSCount;

    /** Field description */
    private final int[] RRCount;

    /** Field description */
    private final int[] RACount;

    /** Field description */
    private final int[] ASCount;

    /** Field description */
    private final int[] ARCount;

    /** Field description */
    private final int[] AACount;

    /** Field description */
    private final int r;

    /** Field description */
    private final int[][][] data;

    /** Field description */
    private final int[][] probeOrders;

    /** Field description */
    private final int[] ch_ids;

    /** Field description */
    private boolean dataLoaded;

    /**
     * Constructs ...
     *
     * @param ch_ids description
     * @param probeOrders description
     */
    public RepDNATests(int[] ch_ids, int[][] probeOrders) {

        if (ch_ids.length != probeOrders.length) {

            System.out.println("from: " + this.getClass().getName() + "constructor");
            System.out.println("# of ch_ids is different from the # of probe orders supplied");
            System.out.print("ch_ids length : " + ch_ids.length + "probe orders: "
                             + probeOrders.length);
            System.exit(0);
        }

        this.ch_ids      = ch_ids;
        this.probeOrders = probeOrders;
        r                = ch_ids.length;
        data             = new int[r][m][m];
        SSCount          = new int[r];
        SRCount          = new int[r];
        SACount          = new int[r];
        RSCount          = new int[r];
        RRCount          = new int[r];
        RACount          = new int[r];
        ASCount          = new int[r];
        ARCount          = new int[r];
        AACount          = new int[r];
    }

    /**
     * Method description
     */
    private void loadData() {

        for (int i = 0; i < r; i++) {

            int ch_id = ch_ids[i];
            Clone[] minimalPathClones = PhysicalMapService.getMinimalPathClones(ch_id,
                                            probeOrders[i], new LinkingCloneStrategy());

            for (int cloneIndex = 1; cloneIndex < minimalPathClones.length; cloneIndex++) {

                Clone clone1   = minimalPathClones[cloneIndex - 1];
                Clone clone2   = minimalPathClones[cloneIndex];
                int clone1Type = clone1.getCloneType();
                int clone2Type = clone2.getCloneType();

                if ((clone1Type == Clone.S_CLONE) && (clone2Type == Clone.S_CLONE)) {
                    SSCount[i]++;
                }

                if ((clone1Type == Clone.S_CLONE) && (clone2Type == Clone.R_CLONE)) {
                    SRCount[i]++;
                }

                if ((clone1Type == Clone.S_CLONE) && (clone2Type == Clone.A_CLONE)) {
                    SACount[i]++;
                }

                if ((clone1Type == Clone.R_CLONE) && (clone2Type == Clone.S_CLONE)) {
                    RSCount[i]++;
                }

                if ((clone1Type == Clone.R_CLONE) && (clone2Type == Clone.R_CLONE)) {
                    RRCount[i]++;
                }

                if ((clone1Type == Clone.R_CLONE) && (clone2Type == Clone.A_CLONE)) {
                    RACount[i]++;
                }

                if ((clone1Type == Clone.A_CLONE) && (clone2Type == Clone.S_CLONE)) {
                    ASCount[i]++;
                }

                if ((clone1Type == Clone.A_CLONE) && (clone2Type == Clone.R_CLONE)) {
                    ARCount[i]++;
                }

                if ((clone1Type == Clone.A_CLONE) && (clone2Type == Clone.A_CLONE)) {
                    AACount[i]++;
                }
            }
        }

        for (int i = 0; i < r; i++) {

            data[i][0][0] = SSCount[i];
            data[i][0][1] = SRCount[i];
            data[i][0][1] = SACount[i];
            data[i][1][0] = RSCount[i];
            data[i][1][1] = RRCount[i];
            data[i][1][1] = RACount[i];
            data[i][2][0] = ASCount[i];
            data[i][2][1] = ARCount[i];
            data[i][2][1] = AACount[i];
        }

        dataLoaded = true;
    }

    public void runIndependenceTest() {

        if (!dataLoaded) {
            loadData();
        }

        IndependenceTest indTest = new IndependenceTest(data);
        double pValue =
            1.0 - new ChiSqrDistribution(indTest.getDF()).cumulative(indTest.getChiSqValue());

        System.out.println("Chi Square: " + indTest.getChiSqValue());
        System.out.println("Chi square DF: " + indTest.getDF());
        System.out.println("P Value: " + pValue);
    }

    public void runHomogeneityTest() {

        if (!dataLoaded) {
            loadData();
        }

        HomogeneityTest homogeneityTest = new HomogeneityTest(data);
        double pValue = 1.0
                        - new ChiSqrDistribution(homogeneityTest.getDF()).cumulative(
                            homogeneityTest.getChiSqValue());

        System.out.println("Chi Square: " + homogeneityTest.getChiSqValue());
        System.out.println("Chi square DF: " + homogeneityTest.getDF());
        System.out.println("P Value: " + pValue);
    }

    /**
     * Class description
     *
     * @version Enter version here..., 12/11/23
     * @author Susanta Tewari
     */
    private class HomogeneityTest {

        double chiSq = -1;


        // transition frequency for each sample (the first dimension)

        /** Field description */
        private int[][][] data;

        /** Field description */
        private double[][] tranMatrix;
        int r;
        int m;

        /**
         * Constructs ...
         *
         * @param data description
         */
        public HomogeneityTest(int[][][] data) {

            if (ArrayUtil.isRagged(data)) throw new Error("data in homogeneity test is ragged");

            this.data  = data;
            r          = data.length;
            m          = data[0].length;
            tranMatrix = new double[m][m];
        }

        public double getChiSqValue() {

            if (chiSq != -1) return chiSq;

            int[][] sampleRowSums = new int[r][m];

            for (int h = 0; h < r; h++) {

                sampleRowSums[h] = FrequencyUtil.getRowSums(data[h]);
            }

            int[] rowSums = SampleTotal.getSampleTotal(sampleRowSums);


            // by summing accross the sample dimension
            int[][] compressedSample = SampleTotal.getSampleTotal(data);

            for (int i = 0; i < m; i++) {

                for (int j = 0; j < m; j++) {

                    tranMatrix[i][j] = compressedSample[i][j] / (double) rowSums[i];
                }
            }

            for (int h = 0; h < r; h++) {

                for (int i = 0; i < m; i++) {

                    for (int j = 0; j < m; j++) {

                        double exp = sampleRowSums[h][i] * tranMatrix[i][j];

                        if (exp > 0) chiSq += Math.pow(data[h][i][j] - exp, 2.0) / exp;
                    }
                }
            }

            if (debug) {

                System.out.println("Sample Row sums");
                LogUtil.print(sampleRowSums);
                System.out.println("Transition Matrix");
                LogUtil.print(tranMatrix, 2, 1, new DecimalFormat("0.####"));
            }

            return chiSq;
        }

        public int getDF() {
            return r * (m - 1) * m;
        }

        public double[][] getTranMatrix() {
            return tranMatrix;
        }
    }

    /**
     * Class description
     *
     * @version Enter version here..., 12/11/23
     * @author Susanta Tewari
     */
    private class IndependenceTest {

        double chiSq = -1;


        // compressed sample or single sample transition table

        /** Field description */
        private int[][] data;
        int m;

        /**
         * Constructs ...
         *
         * @param data description
         */
        public IndependenceTest(int[][][] data) {

            if (ArrayUtil.isRagged(data)) throw new Error("data in homogeneity test is ragged");

            this.data = SampleTotal.getSampleTotal(data);
            m         = this.data.length;
        }

        public double getChiSqValue() {

            if (chiSq != -1) return chiSq;

            int[] rowSums = FrequencyUtil.getRowSums(data);
            int[] colSums = FrequencyUtil.getColSums(data);
            int N         = MathUtil.Sum(rowSums);

            for (int i = 0; i < data.length; i++) {

                for (int j = 0; j < data[i].length; j++) {

                    double exp = (rowSums[i] * colSums[j]) / (double) N;

                    if (exp > 0) {
                        chiSq += Math.pow((data[i][j] - exp), 2.0) / exp;
                    }
                }
            }

            if (debug) {

                System.out.println("Row sums");
                LogUtil.print(rowSums);
                System.out.println("Col sums");
                LogUtil.print(colSums);
            }

            return chiSq;
        }

        public int getDF() {
            return (m - 1) * (m - 1);
        }
    }
}
